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Abstract 



In this paper, we study an adaptive finite element method for a class of a nonlinear 
eigenvalue problems that may be of nonconvex energy functional and consider its ap- 
plications to quantum chemistry. We prove the convergence of adaptive finite element 
approximations and present several numerical examples of micro-structure of matter cal- 
culations that support our theory. 
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1 Introduction 

In this paper, we study adaptive finite element approximations for a class of nonlinear eigen- 
value problems: Find A G M and u G Hq{Q) such that 



where ft C M^, Z G N, a e (0,oo), "1/ : 1] ^ R is a given function. A/" maps a nonnegative 
function over Q to some function defined on Q. 

Many physical models for micro-structures of matter are nonlinear eigenvalue problems of 
type (|l.ip , for instance, the Thomas- Fermi- von Weizsacker (TFW) type orbital- free model used 
for electronic structure calculations |T7l [30l HQ] and the Gross-Pitaevskii equation (GPE) de- 
scribing the Bose-Einstein condensates (BEG) [^Hl]- In the context of simulations of electronic 
structure calculations, the basis functions used to discretize (jl.ip are traditionally plane wave 
bases or typically Gaussian approximations of the eigenfunctions of a hydrogen-like operator. 
The former is very well adapted to solid state calculations and the latter is incredibly efficient 
for calculations of molecular systems. However, there are several disadvantages and limitations 
involved in such methods. For example, the boundary condition does not correspond to that of 
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an actual system; extensive global communications in dealing with plane waves reduce the effi- 
ciency of a massive parallelization, which is necessary for complex systems; and the generation 
of a large supercell is needed for non-periodic systems, which certainly increases the computa- 
tional cost. The finite element method uses local piecewise polynomial basis functions, which 
does not involve problems mentioned above and has several advantages. Although it uses more 
degrees of freedom than that of traditional methods, strictly local basis functions produce well 
structured sparse Hamiltonian matrices; arbitrary boundary conditions can be easily incorpo- 
rated; more importantly, since ground state solutions oscillate obviously near the nuclei, it is 
relatively straightforward to implement adaptive refinement techniques for describing regions 
around nuclei or chemical bonds where the electron density varies rapidly, while treating the 
other zones with a coarser description, by which computational accuracy and efficiency can 
be well controlled. Thus it should be natural to apply adaptive finite element methods to 
solve nonlinear eigenvalue problems resulting from modeling electronic structures. Indeed the 
adaptive finite element method is a powerful approach to computing ground state energies and 
densities in quantum chemistry, materials science, molecular biology and nanosciences [5]'8]. 

The basic idea of a standard adaptive finite element method is to repeat the following 
procedure until a certain accuracy is obtained: 

Solve Estimate Mark Refine. 

Adaptive finite element methods have been studied extensively since Babuska and Rheinboldt 
[3] and have been successful in the practice of engineering and scientific computing. In par- 
ticular, Dorfler ^22j presented the first multidimensional convergence result, which has been 
improved and generalized, see, e.g., (B] [9l |31| |32] [33l |34| |37] for linear boundary value problems, 
Ha [HI [211 [Ml [38] for nonlinear boundary value problems, and [IHl [El [23l [24l [2F for linear 
eigenvalue problems. To our best knowledge, there has been no work on the convergence of 
adaptive finite element approximations for nonlinear eigenvalue problems, though some a pri- 
ori error analyses of finite dimensional Galerkin discretizations for such nonlinear eigenvalue 
problems have been shown in [10] [TT | [T5 l [29 l ET | 142 . 

In this paper, we shall present a posteriori error analysis of an adaptive finite element 
method for a class of nonlinear eigenvalue problems and prove that the adaptive finite el- 
ement algorithm will produce a sequence of approximations that converge to exact ground 
state solutions. As an illustration, we shall also report several numerical experiments on elec- 
tronic structure calculations based on the adaptive finite element discretization [5l |8l [17] , which 
support our theory. Since the nonlinear term occurs, especially the nonlocal convolution in- 
tegration part, there are several serious difficulties in the numerical analysis. Moreover, the 
associated energy functional for this type of problems is usually nonconvex, which brings seri- 
ous difficulties. In our analysis, we shall apply some nonlinear functional arguments and special 
techniques to deal with the local and nonlocal terms carefully. 

This paper is organized as follows. In the coming section, we give an overview of the 
nonlinear eigenvalue problem. In Section [S] we describe the finite element discretization and 
give an a posteriori error analysis. In Section [4] we design an adaptive finite element algorithm 
and prove the convergence of the adaptive finite element algorithm. In Section [5] we show 
some numerical results for micro-structure computations to support our theory. Finally, we 
give several concluding remarks. 

2 Preliminaries 

Let ^7 C be a polytypic bounded domain. We shall use the standard notation for Sobolev 
spaces W^'^{Q) and their associated norms and seminorms, see, e.g., pj 118]. For p = 2, we 
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denote H'{n) = W'^'^{n) and H^{ft) = {v e H^^iQ) : v \dn= 0}, where v \dQ= is understood 
in the sense of trace, || • \\s^q = || • ||s,2,ri, and (•, •) is the standard inner product. The space 
H~^{Q)^ the dual of Hq{Q)^ wih also be used. For convenience, the symbol < will be used in 
this paper. The notation A < B means that A < CB for some constant C that is independent 
of mesh parameters. We shall use Pol{p^ (ci, C2)) to denote a class of functions that satisfy the 
growth condition: 

Pol{p^ (ci, C2)) = {/ • 3 ai, a2 G M such that 

ci\t\P^ai<f{t)<C2\t\P^a2 Vt>0} 

with ci G M and C2,p G [0, 00). 

The weak form of (|l.ip reads as follows: Find A G M and u G Hq{Q) such that 



(2.1) 



( a{Vu,Vv) + {Vu + Af{u^)u,v) = X{u,v) W v e H^{n), 
[ Mln = Z. 

For convenience, we divide nonlinear term JV into local and nonlocal parts: 

where p = u^^ Mi : [0, 00) ^ R is a given function dominated by some polynomial, and A/2 is 
represented by a convolution integration 

A/'2(P) = P''-' / p'^{y)K{--y)dy 



for some given function K and g G M. 

The associated energy functional with respect to this nonlinear eigenvalue problem is ex- 
pressed by 

E{u)= [ {a\\/u{x)\^ ^V{x)u\x)^£{u\x)))dx^^DK{u^'',u^''), (2.2) 
Jn 2^ 

where £ : [0, 00) ^ R is associated with A/i: 

£{s) = [ Mi{t)dt, 
Jo 

and Dk{-^ •) is a bilinear form defined by 

DK{f,g)= I I f{x)g{y)K{x -y)dxdy. 



The ground state solution of problem (jl.ip is obtained by minimizing energy functional (j2.2p 
in the admissible class 

In our discussion, we assume that 

(i) V e L^{n). 

(ii) £ e Pol{p, (ci,C2)) with one of the following conditions: 
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1. ci G (0, oo); 

2. [0,4/3]; 

3. ci G (— oo, 0), p G (4/3, oo) and 



\^ZP-' < inf ( [ |Vii|V / 



2p 



(2.3) 



(iii) A/i(t) G Pol(pi, (ci, C2)) for some pi € [0,2) and Af{{t) e -Po?(p2, (ci,C2)) for some 

P2 e [0,1). 

(iv) G L^(l^), where = {x — y : G r^}. Moreover, K is some nonnegative even 
function and q G [1,3/2). 

Note that these assumptions are satisfied by typical physical models for micro-structures of 
the matter (see, e.g., [Zl |8l HJl |20l |30] ) and condition (|2.3p was first appeared in [7]. 

It is known that under assumptions (i)-(iv), there exists a nonnegative minimizer of energy 
functional (j2.2p [l5j [30l |42] . Moreover, E{v) is bounded below over A under these assumptions 
[iSl |42] , namely, there exist positive constants C and b such that 

E{v)>C-'\\\/v\\l^-b, V^gA (2.4) 

In general, however, the uniqueness of the nonnegative ground state solution is unknown, of 
which the main reason is that energy functional (j2.2p is nonconvex with respect io p = v? for 
almost all molecular models of practical interest. As a result, we introduce the set of ground 
state solutions by 

U = {ueA: E{u) = min£;(v)}. (2.5) 

If It G is a ground state solution, then there exists a corresponding Lagrange multiplier A G M 
such that (A, u) solves (j2.ip and satisfies 

Z\ = E{u)^ [ {Xi{u\x))u\x)-S{u^{x)))dx^{l-^)DK{u^'',u^''). 

We define the set of ground state eigenvalues by 

A = {A G M : (A, ix) solves (|2l]), u G U}. 
The following estimate of the nonlinear term will be used in our analysis. 

Lemma 2.1. Letx^w G H^{Q) satisfy \\x\\i,n-\- ^ C for some constant C. Then there 

exists a constant C > depending on C such that 

imx')x-mw'')w)v<C\\x-w\U,n\\v\\i,i, ^veH^iQ). (2.6) 

Proof We first prove that (|2.6p holds when JV is replaced by local term A/i. Since there exists 
6 e [0,1] such that 

M(x')x - Mr{w')w = {M,{e) + 2eMi{e)){x - ^) 

with ^ = X + ^('^ ~ x)^ we have 

f {M,{x^)x-Mi{w^)w)v= [ (M(a + 2m(e'))(x-H)« ^^veHli^). 
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From assumption (iii), we have that for ah v G HliQ)^ there hold 

I me){X-w)v < ||?'^M|0,3M,qIIX-"'I|0,6/(6-2p.),"IIH|0,6,Q 

Jn 

< m\l%n\\x-M\iMi,a 

and 

/ 2eM[{e){x-w)v < \\e'"^'\\o,3/2,n\\x-w\\o,e,n\\v\\o,6,n 
Jq 

< IICIIo!'6^'llx-«'lli.olkl|i.n, 
where the Holder inequality and the Sobolev inequality are used. Note that 

U\\o,6,a < \\x\\o,6,n + lkl|o,6,o ^ llxlli.n + Iklli.n < C, 

SO we get 

/ mx')x-f^dw')w)v < {ml%n + mlTn)\\x-M\iMi,n 

JQ 

< \\x-wh,n\\v\\i,n V«Gffo'(^)- (2-7) 

For nonlocal term A/2, we obtain from assumption (iv), the Young's inequality, and the 
Holder inequality that 

\\K*{x''-w'^)\\o,oo,n < \\K\\,^dx''-w'no,n 

< ll^llo,Qlle'''-io.6/5.n||x-«^l|o.6.o 

< llx-w'lli.n- 

Hence, for all v G ifQ(il) we have 

< IIX - w||l,Q||xllo'2(2g-l),Qll^llo." 

< \\x-wh,n\\vh,n, (2.8) 
where (7 G [1, 3/2) as in assumption (iv). Similarly, for all v G Hq{Q)^ there holds 

Jn 

< ||^^^"^||o,3/(g-l),^7||X - Uj\\o,6/{5-2q),Q\\v\\l,n 

< \\x-w\\i,n\\v\\i,n. (2.9) 

Taking 1^^, ([221), and identity 

A/'2(x')x - ^2{w^)w = K * (x'^ - ^'^)x''"' + K * ^'^(x'""' - ^''"') 
into account, we obtain 

{mx^)x-mw^)w)v<\\x-w\\iM\i,n ^veHliO), 
which together with (|2.7p leads to (|2.6p . This completes the proof. □ 
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3 Finite element discretizations 

Let {Th} be a shape regular family of nested conforming meshes over Q with size h: there 
exists a constant 7* such that 

Pt 

where, for each T e Th^ hr is the diameter of T, pt is the diameter of the biggest ball contained 
in T, and h = max{/iT : T e Th}- Let £h denote the set of interior faces of Th- And we shall 
also use a slightly abused notation that h denotes the mesh size function defined by 

h{x) = hT, X eT y T eTh. 

Let S^{Q) C H^{ft) be a corresponding family of nested finite element spaces consisting of 
continuous piecewise polynomials over Th of fixed degree n > 1 and Sq{Q) = S^{Q) D i^o(l]). 

Set Vh = sl^{n)nA. 

Under assumptions (i)-(iv), we can obtain the existence of nonnegative ground state solu- 
tions in Vh (see, e.g., We do not have any uniqueness result for this discrete problem 
since the energy functional and the admissible set are nonconvex. We define the set of ground 
state solutions in Vh by 

Uh = {uh e Vh : E{uh) = min E{v)}. (3.10) 

veVh 

We have from (|2.4p that is uniformly bounded: 

sup \\uh\\i,n<C (3.11) 

h<l,Uh€Uh 

with some constant C. 

It is seen that a minimizer Uh G Uh solves 



(3.12) 



with the corresponding finite element eigenvalue A/,, G M satisfying 

Z\h = E{uh) + / {Mi{ul{x))ul{x) - 8{ul{x))) dx^{l- ^)Dk{uI\uI^). (3.13) 

Define 

A/^ = {A/, G R : {Xh.Uh) solves ^^.Uh e Uh}. 

A priori error analysis for (|3.12p has been shown in ^15^. To carry out a posteriori error 
analysis, we need the following result. 

Lemma 3.1. There hold 

hTW{ul)uh\\o,T<\\uh\\i,T yreTh- (3.14) 

Proof. It is obvious that 



hT\mul)Uh\\o,T < \\Uh\\ 



1,T 
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holds for pi = in assumption (iii). By the Holder inequality and the inverse inequality, 

hTW,{ul)Uhh,T < /lT||M(t*|)||0,3,r||w/»||0,6,r < hT\W\\l%^^T\Whfi,T 
< h^T~'''\\u,,\\l';^\\Uh\\l,T<\\UH\\l,T 

for pi G [1, 2) and 

hT\\J\fl{ul)Uh\\o,T < /^T||^i^^'||o,3M,T||^/i||l,6/(3-2pi),T < hh\\l,T 

for pi G (0, 1). Combining with the estimate of A/2 as follows 

hT\\J^2iul)Uhh,T < hT\\K*ul^\\o,^,T\\ur'h,T 

< hT\\K\\,Jumo,n\\u4l'T' < hhlkr, 

where assumption (iv) is used, we obtain (|3.14p . This completes the proof. □ 

Let T denote the class of all conforming refinements by the bisection of an initial triangu- 
lation To- For 7/i G T and any G Uh we define element residual 7^t('^/i) and jump residual 
Je{uh) by 

^T{uh) = XhUh + aAuh - Vuh - J\f{ul)uh in T G 7^, 

Je{uh) = a\/uh\Ti • nt + a\/uh\T2 • ^ = [[a\/uh]]e • nt on e e £h, 

where Ti and T2 are elements in Th which share e and is the outward normal vector of Ti 
on £' for z = 1, 2. Let 0Jh{e) be the union of elements which share e and 0Jh{T) be the union of 
elements sharing a side with T. 

For T eTh, define local error indicator r]h{uh^T) by 

ril{uh,T) = hUnT{uh)\\lT + E he\\Je{uh)\\l,. (3.15) 

eeSh,eCdT 

Given a subset C 1^, we define error estimator r]h{uh^^) by 

Terh,TCuj 

The following result will be used in our convergence analysis though it looks rough. 
Proposition 3.1. Let 7/i G T. If (Xh^Uh) is a solution of then 

rih{uh,T)<\\uh\\i,^^iT) '^T€% 

and 

where the uniform constant > depends only on the data and the mesh regularity. 
Proof We first analyze the element residual. Note that 

hT\\1ZT{uh)\\o,T = hrWXhUh ^ (^Auh - Vuh - M{ul)Uh\\Q^T 

< hrWuhWo^T + hT\\Auh\\o,T + hT\\Vuh\\o,T + hT\\J\f{ul)uh\\o,T^ 
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Using the inverse inequality, assumption (i) and (|3.14p , we have 

hT\\UT{Uh)\\o,T < \\Uh\\l,T^ 

to which similar estimates are true when T is replaced by any T' G uJh{T). 
For the jump residual, by the definition of Je{uh) and the trace inequality, 

hl^'^\\Je{uh)\\o,e = hl^'^\\aVuh\j.^ ■ ni ^ aVuh\j.^ • n^Uce 

< hy\\\\/uM0,e^\\^Uh\TA0,e) 

< hh\\l,UJHiT)' 

Hence 

J2 hl/'\\JeM\\o,e<\\Uh\\l,^,^Ty 

From ()3.1ip . the definition of r]h{uh^T) and rfhiuh-,^)^ we get the desired result. □ 

To present upper and lower error bounds, we introduce an oscillation osCh{uh^T) for any 
T G by 

osc^K,T) = 4||7^TK)-7^TK)||o,T, 

where IZriuh) G Pn-i denotes the projection of IZriuh)- For a subset a; C 11, we define 

TeTh,TCu; 

We have have a standard argument that (see Appendix for a proof) 

Theorem 3.1. Let {X^u) be a regular ground state solution of (|2.ip . If {Xh^Uh) is sufficiently 
close to (A, u), then 

r]h{uh,^) - osch{uh,^) < |A - Xh\ + \\u - Uh\\i,n < r]h{uh,^) ^osch{uh,^). 

Remark 3.1. This result provides the standard upper and lower bounds of the error with 
respect to the error estimator. However^ the hypothesis that (X^u) is a regular solution is 
somehow strong, which can not be proved for most of the problems of practical interest (c.f., 
e.g., Appendix). Anyway, it will not be used in our convergence analysis. 

We define global residual R^(ii^) G H~^{Q) as follows 
(R/,K), v) = Xh{uh, v) - {aVuh, Vv) - {Vun, v) - {J\f{ul)uh, v) \f v e H^{n) 
and see that 

The global residual can be estimated by the local error indicators in the following sense. 
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Theorem 3.2. // {\h,Uh) e^xVh is a solution of \3.12\) . then 

TeTh 

Proof. Let v G Hl{Q) and Vh G Sq{Q) be the Clement interpolant of v satisfying 

11^ - ^/i||o,T < ^T||Vv||o,a;^(T) and \\V - VhW^^ST < ^ ll^^llo,a;^(T) • 

Due to {Ilh{uh)^Vh) = 0, we obtain 

\{Rh{uh),v)\ = \{Rh{uh),v - Vh)\ 

< Yl II^tK)I|0,t||^ - Vhh,T + Yl ll^eK)||0,e||^ - Vh\\o,e 
TeTh \ eeSh,eCdT 



< 



< 



Y II^^TK)||0,T|b||l,..(T)+ E ll^'^'^eK)||0,e||^||l,..(e) 
TeTh \ eeSh,eCdT 

TeTh 



This completes the proof. □ 

4 Convergence of adaptive finite element computations 

We shall first recall the adaptive finite element algorithm. For convenience, we shall replace 
the subscript h (or h^) by an iteration counter k of the adaptive algorithm afterwards. Given 
an initial triangulation To, we can generate a sequence of nested conforming triangulations Tk 
using the following loop: 

Solve Estimate Mark Refine. 

More precisely, to get 7fc+i from Tk we first solve the discrete equation to get Uk on Tk- The 
error is estimated by any Uk G Uk and used to mark a set of elements that are to be refined. 
Elements are refined in such a way that the triangulation is still shape regular and conforming. 

Here, we shall not discuss the step "Solve", which deserves a separate investigation. We 
assume that solutions of finite dimensional problems can be solved to any accuracy efficiently. 
The procedure "Estimate" determines the element indicators for all elements T ^ Tk- A 
posteriori error estimators are an essential part of this step, which have been investigated in 
the previous section. In the following discussion, we use rjkiuk^T) defined by ()3.15p as the 
a posteriori error estimator. Depending on the relative size of the element indicators, these 
quantities are later used by the procedure "Mark" to mark elements in Tk and thereby create 
a subset of elements to be refined. The only requirement we make on this step is that the set 
of marked elements M.k contains at least one element of Tk holding the largest value estimator 
[23l [24] . Namely, there exists one element T^^^ ^ M.k such that 

7?feK,T^^"^) = max7^feK,T). (4.1) 

It is easy to check that the most commonly used marking strategies, e.g.. Maximum strategy, 
Equidistribution strategy, and Dorfier's strategy fulfill this condition. Finally, the marked 
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elements are refined to force the error reduction by ttie procedure "Refine" . The basic algorithm 
in this step is the tetrahedral bisection, with the data structure named marked tetrahedron, 
the tetrahedra are classified into 5 types and the selection of refinement edges depends only 
on the type and the ordering of vertices for the tetrahedra [2 . Note that a few more elements 
T G Tfc \ M-k are partitioned to maintain mesh conformity. It is worth mentioning that we do 
not assume to enforce the so-called interior node property. 

The adaptive finite element algorithm without oscillation marking is stated as follows: 

Algorithm 4.1. 

1. Pick any initial mesh To, and let k = 0. 

2. Solve the system on Tk to get discrete solutions Uk- 

3. Choose any Uk G Uk and compute local error indictors r]k{uk^T) \/ T eTk- 
4- Construct Mk cTk by a marking strategy that satisfies 

5. Refine Tk to get a new conforming mesh Tk+i- 

6. Let k = k and go to 2. 

The purpose of this paper is to prove that Algorithm 14. II generates a sequence of adaptive 
finite element solutions which converge to some ground state solutions of (|2.5p . More precisely, 
we shall prove that 

lim di^iH^iUkM) = 0^ 

k^oo 

lim dist(A/e, A) = 0, 

k^oo 

where 

distj:/! (F, G) = sup inf 

for any F,G C H^{n), and 

dist(74, B) = sup inf \a — b\ 

aeAbeB 

for any A.BcR. 

We first show that the adaptive finite element approximations are convergent. Given an 
initial mesh To, Algorithm 14 . 1 1 generates a sequence of meshes 7i, 72, • • • , and associated discrete 
subspaces 



where 5*00 (^) = I^Sq''{Q) . It is obvious that 6*00 (^) is a Hilbert space with the inner 
product inherited from Hq{Q)^ and there holds 

lim inf \\vk - Voo\\i,Q = V (Q). (4.2) 

We set Ko = Soo{^) n A. 
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Under assumptions (i)-(iv), the existence of minimizers of energy functional (|2.2p in Vex 
can be obtained. Similar to ()2.5p and ()3.1Qp . we introduce the set of minimizers by 

Uoc = {ueVoo: E(u) = min E{v)}. 

veVoo 

We see that Uoo G Uoo solves 



with the corresponding eigenvalue Aoo ^ ^ satisfying 



(4.3) 



ZAoo = ^(^ioo) + / {N'i{uUx))uUx)-S{uUx)))dx^{l-^)^^^ (4.4) 
and we define 

Aoo {Aoo ^ ^ ' (Aoo,^oo) solves Uoo ^ Uoo}' 

Theorem 4.1. If {Uk}ken is the sequence of adaptive finite element approximations generated 
by Algorithm \4.1\ then 

lim Ek — min E{y)^ 
lim d\^iH^{UkMoo) = 0, 
where Ek = E{v){v G Z///c). Moreover, there holds 

lim dist(A/e, Aoo) = 0- 

k^oo 

Proof Following [4T1|42] (see also [15]), let ^i/c G Uk be such that (A/c, ^x/c) solves (|3.12|) in R x 14 
for /c = 1,2,---, and {i^/cmlmeN be any subsequence of {iife}/cGN with 1 < /ci < /c2 < • • • < 

km < • • • . 

Note that (|3.1ip and the Banach-Alaoglu Theorem yield that there exist a weakly conver- 
gent subsequence {uk^,}j^n and i^oo ^ 'S'oo(^) satisfying 



t%rr 



Uoo in H^ift), (4.5) 



we need only to prove 



and 



J^OO 

where (A/e^ . , Uk^ . ) solves (|3.12p and (Aoo, '^^oo) solves (|4.3p . 



E{uoo) = min £^(^), (4.6) 



lim ||ix/c^ -^ioolli,^ = 0, (4.7) 



lim |Afe„. - Aoo I =0, (4.8) 



11 



Since Hq{Q) is compactly imbedded in L^{Q) for p G [2,6), by passing to a further subse- 
quence, we may assume that Uk^. Uoo strongly in Lp{Q) as j ^ oo. Thus we can derive 

lim [ £{ul^ )= [ £{ul), 
lim Dk{uII ^,ull) = Dk{u^^, u^£). 

and hence 

liminf^K^.) >^(2ioo). (4.9) 

Note that (|4.2p implies that {uk^ ,} is a minimizing sequence for the energy functional in 
5*00 (^), which together with (|4.9p and the fact that {uk^.} converge to Uoo strongly in L'^{Q) 
leads to Uoo G Uoo^ namely, 

lim E{uk^,) E{uoo) min E{v). 

j^OO ^ veVoo 

Consequently, we obtain that each term of E{v) converges and in particular 

lim \\VUk^,\\o,Q = \\VUoo\\o,Q' 
j^oo ^ 

Using (|4.5p and the fact that Hq{Q) is a Hilbert space under the norm ||V • \\o,n^ we have 

lim II -Uoc)\\o,Q = 0, 

which implies (|4.7p 

Using (I3.13p . (|44|) . (|46|) and (j4?7j), we immediately obtain (j4]8j). This completes the proof. 

□ 

Following the ideas in [23l [24l [34l |36] , we then prove the convergence of the a posteriori 
error estimators and the weak convergence of residual R/c(iXfe), which will be used to prove that 
the adaptive finite element approximations converge to the ground state solutions. Given the 
sequence {Tk}keN^ for each k eN we define 

= {T eTk-.T eTu^ l>k} and 7^ = Tk\ r+. 

Namely, Tj^ is the set of elements of Tk that are not refined and Tj^ consists of those elements 
which will eventually be refined. Set 

= ^Ter+^k{T) and ft^ = Urer^^kiT). 

Note that the mesh size function hk = hkix) associated with Tk is monotonically decreasing 
and bounded from below by 0, we have that 

^oo(^) lim hk{x) 

k^oo 

is well-defined for almost all x G 1^ and hence defines a function in L^(yt). Moreover, the 
convergence is uniform |34] . 
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Lemma 4.1. If {hk}keN is the sequence of mesh size functions generated by Algorithm \4 1\ 
then 

lim \\hk - /ioo||o,oo,^7 

and 

lim \\hkXnlh.oo,n = 0, 
where Xqo is the characteristic function of 

Lemma 4.2. If {uk}ken is a sequence of ground state solutions of \3.12]) in {Vk}ken obtained 
by Algorithm \4-.l\ then 

lim max rikiuk.T) = 0. 

k^ooTeMk 

Proof We see from the proof of Theorem 14.11 that for any subsequence {uk^} of {uk}^ there 
exist a convergent subsequence {uk^. } and Uoo G Uoo such that 

Now it is only necessary for us to prove that 

1™ ^i^^x r]k^ {uk , T) = 0. 

In order not to clutter the notation, we shall denote by {ii/c}/cGN the subsequence {uk^.}men^ 
and by {Tk}keN the sequence {7^^^. }mGN- 
Let T/e e M-k be such that 

r]k{uk,Tk) = max r]k{uk,T). 
Using Proposition 13. 1[ we obtain 

r]k{Uk,Tk) < \\Uk\\l^uJk{Tk) < ll^/c-^oo||l,^7+||^oo||l,a;fe(Tfe)- (4.10) 

Since TkeMkC T^, we have 

\^k{Tk)\ < h^T, < \\hkXnl\\loo,Q ^ ^ as k ^ oo, 

where Lemma 14.11 is used. From Theorem 14. If we have that the first term in the right hand 
side of (j4.1Qp tends to zero, too. This completes the proof. □ 

Lemma 4.3. If {uk}keN is a sequence of ground state solutions of \3.12\) in {Vk}ken obtained 
by Algorithm \4-l\ then 

lim {Rk{uk),v) =0 y V e H^{ft). 

k^oo 

Proof Using similar arguments as that in the proof of Theorem l4.H for any subsequence {uk^} 
of {uk}^ there exist a convergence subsequence {uk^. } and Uoo G Uoo such that 
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and we need only to prove 

lim (R/e^.(iife^.),v) =0 y V e H^{n). 



J^OO 

t2( 



Since Hq{Q) is dense in Hq{Q)^ it is sufficient to prove 

lim(Rfe K ),«)=0 \fveHl{a). (4.11) 

For simplicity of notation, we denote by {i^fe}feGN the subsequence {uk^}meN^ and by {Tk}keN 
the sequence {Tk^jmeN- 

Let Vk G Vk be the Lagrange's interpolation of v. Since 

{Rk{uk),Vk) = 0, 

we have from Theorem 13.21 that 



\{Rk{uk),v)\ = \{Rk{uk),v - Vk)\ < ^ r]k{uk,T)\\v - VkWi^ojkiT)- 

TeTk 

Let n G N and k > n. By definition, C Tj^ C Tk- Thus we have 

\{Rk{Uk),v)\ < ^ Vk{Uk,T)\\v - VkWl^uJkiT) ^ Vk{Uk,T)\\v-Vk\\l^oJkiT) 

Ter+ Terk\r;t 

< Vk{Uk,T^)\\v - Vfe||i,n+(T) + Vk{Uk,Tk \ T;^)\\V - V/c||l,QO(T). 

Using Proposition 13.11 we get 

r]k{uk,Tk \ Tn) < r]k{uk,Tk) < Cr^, 
which together with the interpolation estimate yields 

|(RfcK),^)| < (7?fcK,r+)+C^||/lnX^7o||o,oo,Q)|b||2,Q. (4.12) 

Now we shall use (|4.12p to prove (|4.1ip . Let £ > be arbitrary. Lemma 14.11 implies that 
there exists n G N such that 

Cryll/^nXQO ||0,oo,^7 < S. (4.13) 

Since C C Tk and the marking strategy ()4.ip is reasonable, we arrive at 

7?fc(iife,r+) < (#r+)^/' max 7?fc(iife,T) < (#r+)^/' max r]k{uk,T). 

Ter+ TeMk 

By Lemma |4]2l we can select N > n such that 

r]k{uk,T^)<e yk>N. (4.14) 

Thus we obtain (|4TT]) by combining (|4?T2)) . (|4?T3l) and (|4?T4|) . This completes the proof. □ 

Finally, we prove the main result of this paper. 
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Theorem 4.2. Given a sufficiently fine initial mesh To- If {Uk}keN is the sequence of adaptive 
finite element approximations generated by Algorithm \4.1\ then 

lim Ek = mmE{v), (4.15) 
lim distHi{Uk,U) = 0, (4.16) 

k^oo 

and 

lim dist(A/e,A) = 0. (4.17) 

k^oo 

Proof We shall use a similar argument as that in the proof of Theorem 14.11 Let be 
such that {Xk^Uk) solves (j3.12p for k = 1,2,- ••. For any subsequence {uk^} of {uk}^ there 
exist a convergent subsequence {uk^ .} and Uoq G Uoo such that 

Ukm^ ^oo in i^o(^), 

^krrij ~^ -^00 5 

where (Aoc'^^oo) solves (|3.12p . It is only necessary for us to prove that Uoq G which derives 
(j4.15p and (j4.16p directly, and implies (|4.17p by noting (|2.ip and (|3.12p . For simplicity, we 
denote by {ii/c}feGN the convergent subsequence {uk^.}meni and by {Tk}ken the subsequence 

We first prove that limiting eigenpair (Aoc'^^oo) is also an eigenpair of (|2.ip . Note that 

Aoo(^oo, ^) - Q^(Viioo, Vv) - {VUoo + M{ul^)Uoo, V) - (R/e(iife), v) 

= (Aoo^oo - Afeix/c, v) - a{V{uoo - Uk), Vv) 

-(y(uoo - Uk), v) - {JV{ul^)uoo - M(ul)uk, v) y ve H^i^), 

we obtain from (|2.6p that 

|Aoo(^oo,^) - a(Viioo, Vv) - {VUoo ^^f{ul^)Uoo,v) - (Rfc (iifc) , | 

< ||V(^Xoo - Uk)\\o,n\\^v\\o,n + ||^||o,j7||^oo - ^/c||o,3,q||^||o,6,q 

+ 11^00 - Uk\\i,Q\\v\\i^n + {\\uoo - Uk\\o,n + \^k - ^oc\)\\v\\o,n V v G H^{n). (4.18) 

Since Xk Aoo and Uk Uoo in Hq{Q)^ the right hand side of (|4.18p tends to zero when k 
tends to infinity. Using Lemma H]3] and identity 

Aoo(^oo,^) - a(Vlioo, Vv) - (VUoo ^ M{u^^)Uoo,v) 

= Aoo(^oo,^) - a(Viioo, Vv) - {Vuoo M{ul^)uoo,v) - {Rk{uk),v) + {Rk{uk),v), 
we arrive at 

a(Viioo, Vv) + (Vuoo ^ M{uIq)uoo,v) = Xoo{uoo,v) y V e Hq{Q). 

Now we prove that for a sufficiently fine initial mesh, the limiting eigenfunction Uqq is a 
ground state solution. Set 

W = {w e HQ{fl) : is an eigenfunction of (|2.ip }. 
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Note that U C W, the ground state solutions in U minimize energy functional (|2.2p , which is 
continuous over HliQ), we can choose a mesh To such that 

Eo = E(v)< min E(w) W v eUo, 
wew\u 

where the fact 

lim inf \\v-w\\iQ = WweHl,(n) 

is used. Due to 7^ C To, we have Ek < Eq and obtain Uqo G U. This completes the proof. □ 
If we make a further assumption that 

S'\t)>0 fortG[0,oo), (4.19) 

then energy functional 

E{Vp)= I {a\V^p\^^V{x)p{x)^8{p{x)))dx^^DK{p\p'') 
Jn ^Q. 

is strictly convex on convex set {p > : ^/p G A} and hence there exists a unique minimizer 
of ()2.2p in admissible class A. Note that the minimizer of ()2.2p in Vk is unique when initial 
mesh To is fine enough (c.f., e.g., [42 ), we have 

Corollary 4.1. Assume that the hypothesis of Theorem \4.^ and (|4.19p are satisfied. If{X,u) G 
R X ^ i5 the ground state solution of (|2.ip and {Xk^Uk) G M x 14 ^-^ the ground state solution 
of (j3:T2l) . then 

lim \\uk - u\\ 0, 

k^oo 

lim IIAfe - All = 0. 

k^oo 

5 Numerical examples 

In this section, we will report on some numerical experiments for both linear finite elements 
and quadratic finite elements in three dimensions to illustrate the convergence of adaptive finite 
element approximations. 

Our numerical computations are carried out on LSSC-II in the State Key Laboratory of 
Scientific and Engineering Computing, Chinese Academy of Sciences, and our codes are based 
on the toolbox PHG of the laboratory. All of the computational results are given in atomic 
unit (a.u.). 

Example 1. Consider the ground state solution of CPE for BEC with a harmonic oscillator 
potential 

V{x, y, z) = ]plx^ + -ily^ + 7'^'), 

where = 1,7?/ = 2,7^ = 4. We solve the following nonlinear problem: Find (X^u) G 
R X iTo(^) such that ||ii||o,Q = 1 and 

J (-|A + V + /3|iip)2i = Xu in^^, 
1 = on dVt, 
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where (3 = 200 and Q = [-8,8] x [-6,6] x [-4,4]. 

The convergence of energies and the reduction of the a posteriori error estimators are 
presented in Figure [STTl which support our theory and that the a posteriori error estimators are 
efficient. Some cross-sections of the adaptively refined meshes constructed by the a posteriori 
error indicators are displayed in Figure 15.21 




number of degrees of freedom 



number of degrees of freedom 



Figure 5.1: Left: Convergence curves of energy for BEC. Right: Reduction of the a posteriori 
error estimators using linear and quadratic elements. 




Figure 5.2: The cross-sections on z = of adaptive meshes using linear (left) and quadratic 
(right) elements. 

In the next two examples, we shall carry out the ground state energy calculations of atomic 
and molecular systems based on TFW type orbital- free models. The nonlinear term is given 

by 

.2/ 



Jm3 \--y\ 3 



where Ctf = Joi^'^'^)^ ^xc is the exchange-correction potential. The exchange-correction 
potential used in our computation is chosen as 

vM = y^''^{p) + yc''^{ph (5.1) 

where 



1/3 
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0.0311 In r, - 0.0584 + O.OOlSr, In r, - 0.0084r, if r, < 1 

-(0.1423 + 0.0633r, + 0.1748v^)/(l + 1.0529v^ + 0.3334r,)2 if r, > 1 

1/3 



and To 



(-) 

Example 2. Consider the TFW type orbital- free model for helium atoms. The external 
electrostatic potential is V{x) = —2/\x\. Then we have the following nonlinear problem: Find 
{X,u) G M X H^{Q) such that \\u\\l ^ = 2 and 



1 2 f \u(y)\'^ 5 

■--Au- —u^u ^^\ dy^ -Ctfu^^^ ^v^c{u'^)u = Xu inn, 
10 1^1 Jq \x — y\ 6 

u = on dfl, 



where Q = (-5.0,5.0)^. 

The convergence of energies and the reduction of the a posteriori error estimators are 
shown in Figure 15. 3[ which support our theory. The cross-sections of the adaptive meshes are 
displayed in Figure [5l4| from which we observe that more refined meshes (nodes) appear in the 
area where the nuclei are located. 




Figure 5.3: Left: Convergence curves of energy for the helium atom. Right: Reduction of the 
a posteriori error estimators using linear and quadratic elements. 



Example 3. Finally, we consider an aluminum cluster in the face centered cubic lattice 
consisting of 3 x 3 x 3 unit cells with 172 aluminium atoms, where the GHN pseudopotential 
[26] is used. We solve the following nonlinear problem: Find (X^u) G M x Hq^Q) such that 
\\u\\l^ = 172 and 

/ -Lau^V^^.^^'^u^u [ ^-^^ = Xu inn, 

^ 10 ^ Jn\x-y\ 3 

[ u = on dQ, 

where Q = (-25.0,25.0)^. 

The convergence of energies and the reduction of a posteriori error estimators are shown in 
Figure 15.51 The cross-sections of the adaptive meshes are displayed in Figure 15.41 We observe 
that with the a posteriori error estimators, the refinement is carried out automatically at the 
regions where the computed functions vary rapidly, especially near the nuclei. As a result, 
the computational accuracy can be controlled efficiently and the computational cost is reduced 
significantly. 



18 



Figure 5.4: The cross-sections on 2; = of adaptive meshes using hnear (left) and quadratic 
(right) elements. 




Figure 5.5: Left: Convergence curves of energy for the aluminium cluster in FCC lattice. Right: 
Reduction of the a posteriori error estimators using linear and quadratic elements. 

6 Concluding remarks 

We have analyzed adaptive finite element approximations for ground state solutions of a class of 
nonlinear eigenvalue problems. We have proved that the adaptive finite element loop produces 
a sequence of approximations that converge to the set of exact ground state solutions. This 
result covers many mathematical models of practical interest, for instance, the Bose-Einstein 
condensation, the TFW model in the orbital-free density functional theory, and Schrodinger- 
Newton equations in the quantum state reduction [151 EZl [35] where the integration kernel K 
is negative. We have also applied adaptive finite element discretizations to micro- structure 
of matter calculations, which support our theory. It is shown by Figure 15.11 Figure 15. 3| and 
Figure [531 that we may have some convergence rates of adaptive finite element approximations. 
Indeed, it is our on-going work to study the optimal complexity of adaptive finite element 
approximations for such nonlinear eigenvalue problems, which requires a new technical tool 
and will be addressed elsewhere [16 . 

Acknowledgements. The authors would like to thank Dr. Xiaoying Dai, Prof. Lihua 
Shen, and Dr. Dier Zhang for their stimulating discussions and fruitful cooperations on elec- 
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Figure 5.6: The cross-sections on z = of adaptive meshes using hnear and quadratic finite 
elements. 



tronic structure computations that have motivated this work. The authors are grateful to Prof. 
Linbo Zhang and Mr. Hui Liu for their assistance on numerical computations. 



Appendix 

We may follow the framework in [39] to derive an upper bound and a lower bound of the a 
posteriori error estimate. Let 

X = Y = Hl(^l), 

{F{\u),{ii,v)) = / {aVuVv^Vuv^J\f(u^)uv-\uv) ^ ii{ j - Z). 

Jn Jn 

It is seen that (|2.ip may be written as 

F{X,u) = (A.l) 

and F G C^(R X X,R X F*), where F* is the dual space of Y. A solution {X,u) of (|AU is 
said to be regular if the following equation 

DF{\,u)-{fi,v) = 

is uniquely solvable in R x X for each {k, f) gRxY*, where DF{X, u) : R x X -> R x F* is 
the Frechet derivative of F at (A, u). It is seen that (A, w) e R x Hq{Q) is a regular solution of 



(TOIl if 

{{E"{u)-\)v,v)y',x>j\\Vv\\1^ ^veHUn)nu^ (A.2) 
is true for some constant 7 > (see, e.g., [29 J, where 

{{E"{u)- \)v,w) = a{Vv,Vw) + {{V+Af{u'^)-\)v,w)+2{Af{ 
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and 

= {w e L^{n) : {u,w) = 0}. 

It has been proved that ()A.2p is satisfied by some special TFW models that are of convex 
functional (see, e.g., fW\ flT]). 
Let 

= Yh = s^{n) 

and define 

{Fh{Xh.Uh)Af^.v)) = {F{\h,Uh),{^i,v)) y {fi,v) eRx s^{n). 

We see that F/^ G C(M x X/i,R x Y^) and is an approximation of F. Obviously, finite element 
eigenvalue problem (j3.12p is equivalent to 

Fh{Xh.Uh) = 

and (Xh^Uh) is an approximate solution of (jA.ip . 

The following proposition in [39, Section 2.1] yields a posteriori error estimates in the 
neighborhood of (A, li) that satisfies (jA.ip . 

Proposition A.l. Let (X^u) be a regular solution of HA.l]} . If DF is the derivative of F 
and DF is Lipschitz continuous at {X^u), then the following estimate holds for all {Xh^Uh) 
sufficiently close to this solution: 

\\F{Xh,Uh)\\MxY* < \Xh - A| + \\uh -u\\x< \\F{Xh,Uh)\\^y,Y*- (A.3) 

It is shown from (|A.3p that ||F(A^, Uh)\\RxY* is a posteriori error estimator. When we apply 
the general approach in [39^ Sections 3.3-3.4] to 

a{x,u^\/u) = a\/u and b{x^u,\/u) = Xu — Vu — J^{u'^)u 

by taking 

a^(x, Uh, Vuh) = aVuh, bh{x, Uh, Vuh) = Xuh - Vuh - M{u\)uh, 

and 

r]T = r]h{uh,T), St = osch{uh,T), 
we then gives a proof of Theorem 13.11 
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